#'
#'divide a large plot into small nonoverlap quadrats grid
#'
#'@usage
#'quadratize(com, nx=NULL,ny=NULL,qside=NULL,iorder="x")
#'quadrat_number(quad)
#'plot(quad)
#'
#'
#'@param com a community object
#'@param quad a quadrat object returned by the quadratize function
#'@param nx number of quadrat in the x axis
#'@param ny number of qudart in the y axis
#'@param qside a vector of the lengths of quadrat, incompatable with nx and ny. it is only used when nx or ny are missing.
#'@param iorder whether order of quadrat index increase by x-axis or y-axis.
#'
#'@details
#'Order of the quadrat index is controled by the iorder parameter, if iorder="x" (default), then quadrat 
#'index increase along the x-axis first; if iorder="y", then the index increase along the y-axis first.
#'
#'@return
#'a quadrat object with a "quadrati" column indicates the indexes of quadrat each individual belongs to.
#'it also has attributes "qsize" and "iorder" represent the size of each quadrat and order of quadrat indexes, respectively
#'
#'@examples
#'data(BCI)
#'com=quadratize(BCI,5,5,iorder="x")
#'com
#' 
#'com2=quadratize(BCI,qside=c(20,20))
#'com2=quadratize(BCI,qside=c(25,25),iorder="y")
#'com2
#'
#'com3=rcommunity()
#'plot(quadratize(com3,5,5,iorder="x"),showIndex=TRUE)
#'plot(quadratize(com3,5,5,iorder="y"),showIndex=TRUE)
#'
#'
#' @export
quadratize=function(com,nx=NULL,ny=NULL,qside=NULL,iorder="x"){

  if(!is.null(qside)){
    if(length(qside)!=2 | any(is.na(qside)))
      stop("length of qside should be 2, and it can't contain any missing value")
  }
  
  plotdim=attr(com,"plotdim")
  
  if(is.null(nx) & !is.null(qside[1])){
    nx=plotdim[1]/qside[1]
  }
  if(is.null(ny) & !is.null(qside[2])){
    ny=plotdim[2]/qside[2]
  }
  
  if(is.null(qside)){
    qside[1]=plotdim[1]/nx
    qside[2]=plotdim[2]/ny 
  }
  
  if(qside[1]!=round(qside[1],0) | nx!=round(nx,0) | ny!=round(ny,0))
    stop("The whole plot can not be divided by integer number of quadrats with such size")
  
  xi_breaks=seq(-1e-5,plotdim[1]+1e-5,length.out=nx+1)
  yi_breaks=seq(-1e-5,plotdim[2]+1e-5,length.out=ny+1)
  xi=as.numeric(cut(com$x,breaks=xi_breaks))
  yi=as.numeric(cut(com$y,breaks=yi_breaks))
  if(iorder=="x"){
    ploti=xi+(yi-1)*nx
  }else{
    ploti=yi+(xi-1)*ny
  }
  
  com$quadrati=ploti
  attr(com,"qsize")=qside
  attr(com,"iorder")=iorder
  class(com)=c("quadrat",class(com))
  return(com)
}

#number of quadrats
quadrat_number=function(quad){
  plotdim=attr(quad,"plotdim")
  qsize=attr(quad,"qsize")
  return(prod(plotdim/qsize))
}

print.quadrat=function(quad){
 print.community(quad)
 nq=quadrat_number(quad)
 qsize=attr(quad,"qsize")
 cat(paste("divided into ",nq," nonoverlap quadrats with size ",qsize[1],"*",qsize[2],"\n",sep=""))
}

plot.quadrat=function(com,showIndex=FALSE,...){
  p=plot.community(com,...)
  qsize=attr(com,"qsize")
  plotdim=attr(com,"plotdim")
  p=p+geom_vline(xintercept=seq(0,plotdim[1],qsize[1]),colour="grey")+
    geom_hline(yintercept=seq(0,plotdim[2],qsize[2]),colour="grey")
  if(showIndex){
    iorder=attr(com,"iorder")
    tx=seq(qsize[1]/2,plotdim[1],qsize[1])
    ty=seq(qsize[2]/2,plotdim[2],qsize[2])
    labels=1:(quadrat_number(com))
    if(iorder=="x"){
      txx=rep(tx,times=length(ty))
      tyy=rep(ty,each=length(tx))
    }else{
      txx=rep(tx,each=length(ty))
      tyy=rep(ty,times=length(tx))
    }
    p=p+annotate("text",x=txx,y=tyy,label=labels,colour="grey")
  }
  return(p)
}